rm(list = ls())
library(pacman)
p_load(lme4,lmerTest ,performance, tidyverse,lubridate, lattice, devtools, VGAM, ggpubr, plotly, knitr)
rm(list = ls())
setwd("F:/STUDY HCA IN TAIWAN/0. ideas for thesis/COVID19 and NPIs/Data analysis/dataset/data_sep")
load(file="processed_data.RData")
# check the Missing data on test.1000
miss.test <- data0 %>% filter(is.na(test.thousand) ==T) %>%
select(iso_code, country, date.date, test.thousand )
miss.test.country <- data0 %>% filter(is.na(test.thousand) ==T) %>%
distinct(iso_code)
## Countries dont have the missing data
no.mist.test <- data0 %>% filter(is.na(test.thousand) ==F) %>%
distinct(iso_code, country) # 31 countries
save(no.mist.test, file = "no.miss.test.RData")
# Check missing of week.case.eff
miss.case <- data0 %>% filter(is.na(total.cases) == T)%>%
select(iso_code, country, date.date, week.case.eff )
miss.case.cnt <- data0 %>% filter(is.na(total.cases) == T)%>%
distinct(iso_code)
setwd("F:/STUDY HCA IN TAIWAN/0. ideas for thesis/COVID19 and NPIs/Data analysis/dataset/data_sep")
vac.week20 <- data0.anal %>% filter(time == 19) %>% select(iso_code, country, total.cases,tt.vac.100, pp.vac.100, ful.vac.100)
write.table(vac.week20, file="vac_week20_plotdata.csv", sep = ",")
save(vac.week20, file = "vac_week20.RData")
# check.data <- data0.anal %>% select(id, iso_code, time, date.date, tt.vac.100, pp.vac.100, vac.100.lag, ful.vac.100, test.thousand, total.cases, week.case.eff)
#
#
# check.data <- data0.ful %>% select(id, iso_code, time, date.date, tt.vac.100, pp.vac.100, ful.vac.100, test.thousand, avg.string, week.case.eff)
#
# check.data <- vac.period %>% select(id, iso_code, date.date, population, tt.vac.100, pp.vac.100, ful.vac.100, total.cases, test.thousand)
#
# check.vac <- data0.anal %>% filter(is.na(vac.100.cate) == T) %>% select(iso_code, id, time, pp.vac.100, vac.100.cate)
# hist(data0.anal$week.case.eff
# )
# check <- data0.anal %>% filter(week.case.eff <0) %>% select(iso_code, id, week.case.eff)
# check <- data0.anal %>% filter(iso_code == "KAZ") %>%
# select(iso_code, date.date, time, total.cases, week.case.eff
# )
# KAZ data wrong
check.time <- data0.anal %>% select(id, week.case.eff, growth.rate.manual, timeStart, time)
## Adding missing grouping variables: `iso_code`
# Group plot
plot.group.case <- ggplot(data0.anal, aes(time+2, week.case.eff)) +
stat_smooth(method = "loess", se = F, size = 2) +
stat_summary(aes(col = country), fun.y = mean, geom = "line", alpha = 0.5) +
xlab("Week") +
ylab(" Average weekly growth rate")+
theme_bw() +
theme(strip.background = element_blank()) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
theme(legend.position = "bottom") +
scale_x_continuous(breaks = seq(2, 25, by = 1))
## Warning: `fun.y` is deprecated. Use `fun` instead.
plot.group.case
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 3 rows containing non-finite values (stat_smooth).
## Warning: Removed 3 rows containing non-finite values (stat_summary).
plot.group.string <- ggplot(data0.anal, aes(time, string.idx)) +
stat_smooth(method = "loess", se = F, size = 2) +
stat_summary(aes(col = country), fun.y = mean, geom = "line", alpha = 0.5) +
xlab("Week") +
ylab(" Stringency index")+
theme_bw() +
theme(strip.background = element_blank()) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
theme(legend.position = "bottom") +
scale_x_continuous(breaks = seq(0, 23, by = 1))
## Warning: `fun.y` is deprecated. Use `fun` instead.
plot.group.string
## `geom_smooth()` using formula 'y ~ x'
cairo_pdf(filename = "F:/STUDY HCA IN TAIWAN/0. ideas for thesis/COVID19 and NPIs/Data analysis/results/plot_sep/fig1_NPI_vac.period_vac.period_Sep27.pdf",
width = 18, height = 8, #inch
onefile = TRUE, family = "Segoe UI")
merge1 <- ggarrange(plot.group.string, plot.group.case, hjust = -1,
ncol = 2, nrow = 1)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 3 rows containing non-finite values (stat_smooth).
## Warning: Removed 3 rows containing non-finite values (stat_summary).
merge1
dev.off()
## png
## 2
ggplotly(plot.group.case)
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 3 rows containing non-finite values (stat_smooth).
## Warning: Removed 3 rows containing non-finite values (stat_summary).
ggplotly(plot.group.string)
## `geom_smooth()` using formula 'y ~ x'
#https://www.guru99.com/r-scatter-plot-ggplot2.html#2
# clean daily grow rate = 0, and NAs
## NAs = 36 obs
## zero growth rate = 58
## Check zero growth rate
data0.anal1 <- data0.anal
check <- data0.anal1 %>% filter(week.case.eff == 0 )
check1 <- data0.anal1 %>% filter(week.case.eff != 0 ) %>% select(id, week.case.eff)
## Adding missing grouping variables: `iso_code`
## Check NAs
summary(data0.anal$week.case.eff)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.000000 0.002658 0.004951 0.010112 0.009787 0.357987 3
#probit transformation
data0.anal <- data0.anal %>% filter(week.case.eff != 0) %>% mutate(probit.case = probitlink( week.case.eff))
qqnorm(data0.anal$probit.case)
hist(data0.anal$probit.case)
boxplot(data0.anal$probit.case)
summary(data0.anal$test.policy.red)
## Warning: Unknown or uninitialised column: `test.policy.red`.
## Length Class Mode
## 0 NULL NULL
data0.anal <- data0.anal %>% mutate(test.policy.red = case_when(
test.policy==0~0,
test.policy==1~0,
test.policy==2~0,
test.policy==3~1
))
data0.anal$test.policy.red <- factor(data0.anal$test.policy.red)
data0.anal$pop_log <- log(data0.anal$population)
hist(data0.anal$pop_log)
data0.anal$pop_denlog <- log(data0.anal$pop_density)
hist(data0.anal$pop_denlog)
data0.anal$median_age.log <- log(data0.anal$median_age)
hist(data0.anal$median_age.log)
data0.anal$gdp.log <- log(data0.anal$gdp)
hist(data0.anal$gdp.log)
data0.anal$incomegroup1 <- factor(data0.anal$incomegroup, levels = c("Low income", "Lower middle income", "Upper middle income" , "High income"))
data0.anal <- data0.anal %>% mutate(incomegroup2 = case_when(
incomegroup1 == "Low income" ~ "LMIC",
incomegroup1 == "Lower middle income" ~ "LMIC",
incomegroup1 == "Upper middle income" ~ "upper_income",
incomegroup1 == "High income" ~ "high_income"
))
data0.anal$incomegroup2 <- factor(data0.anal$incomegroup2, levels = c( "LMIC", "upper_income","high_income"))
#
# Sig at >30
# data0.anal <- data0.anal %>% mutate(vac.bin = case_when(
# pp.vac.100 <=5 ~ 0,
# pp.vac.100 <=10 & pp.vac.100 >5 ~ 1,
# pp.vac.100 <=30 & pp.vac.100 >10 ~ 2,
# pp.vac.100 > 30 ~ 3,
# ))
# data0.anal <- data0.anal %>% mutate(vac.bin = case_when(
# pp.vac.100 <=1 ~ 0,
# pp.vac.100 <=8 & pp.vac.100 >1 ~ 1,
# pp.vac.100 <=16 & pp.vac.100 >8 ~ 2,
# pp.vac.100 > 16 ~ 3,
# ))
data0.anal$pp.vac.1 <- data0.anal$pp.vac.100
data0.anal <- data0.anal %>% mutate(pp.vac.100 = case_when(
pp.vac.1 <1 ~ 0,
pp.vac.1 <5 & pp.vac.1 >=1 ~ 1,
pp.vac.1 <10 & pp.vac.1 >=5 ~ 2,
pp.vac.1 <30 & pp.vac.1 >=10 ~ 3,
pp.vac.1 > 30 ~ 4,
))
data0.anal$pp.vac.100 <- factor(data0.anal$pp.vac.100)
summary(data0.anal$ful.vac.100)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 0.670 4.685 3.185 60.780
data0.anal <- data0.anal %>% mutate(ful.100.cate = case_when(
ful.vac.100 <=10 & ful.vac.100 >=0 ~ 0,
ful.vac.100 <=30 & ful.vac.100 >10~ 1,
ful.vac.100 <100 & ful.vac.100 >30 ~ 2
))
data0.anal$ful.100.cate <- factor(data0.anal$ful.100.cate)
library(performance)
aic <- function(model){
fix2 <- performance(model)
x<- as.character(substitute(model))
fix2 <- fix2 %>%
mutate( npi = x)
}
model.school.cls <- lmer(probit.case ~ school.close.red +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.school.cls)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ school.close.red + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 87.5 122.1 -35.8 71.5 547
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7392 -0.5252 -0.0208 0.4549 7.4866
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2254225 0.47479
## time 0.0008684 0.02947 -0.63
## Residual 0.0477187 0.21845
## Number of obs: 555, groups: iso_code, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.613e+00 9.421e-02 3.133e+01 -27.734 < 2e-16 ***
## school.close.red2 -2.283e-02 4.198e-02 5.399e+02 -0.544 0.58687
## school.close.red3 1.118e-01 3.445e-02 5.414e+02 3.246 0.00124 **
## time -7.892e-04 5.806e-03 2.765e+01 -0.136 0.89286
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sch..2 sch..3
## schl.cls.r2 -0.169
## schl.cls.r3 -0.233 0.549
## time -0.615 -0.035 -0.002
x1 <- aic(model.school.cls)
model.work.cls <- lmer(probit.case ~ work.close.red +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.work.cls)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ work.close.red + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 104.8 139.3 -44.4 88.8 547
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3332 -0.5075 -0.0527 0.5271 7.2350
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2281360 0.47764
## time 0.0009398 0.03066 -0.63
## Residual 0.0491885 0.22178
## Number of obs: 555, groups: iso_code, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.575427 0.096201 33.054626 -26.771 <2e-16 ***
## work.close.red2 0.034585 0.043643 549.693377 0.792 0.428
## work.close.red3 0.025341 0.047098 547.950575 0.538 0.591
## time -0.001756 0.006040 28.153571 -0.291 0.773
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) wrk..2 wrk..3
## wrk.cls.rd2 -0.271
## wrk.cls.rd3 -0.257 0.670
## time -0.594 -0.072 -0.040
x2<- aic(model.work.cls)
model.pubevent <- lmer(probit.case ~ cancel.pubevent.red +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.pubevent)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ cancel.pubevent.red + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 103.3 133.6 -44.7 89.3 548
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3111 -0.5075 -0.0424 0.5259 7.1903
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2274454 0.47691
## time 0.0009387 0.03064 -0.64
## Residual 0.0492832 0.22200
## Number of obs: 555, groups: iso_code, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.565247 0.098900 36.761623 -25.938 <2e-16 ***
## cancel.pubevent.red2 0.012776 0.043990 546.532929 0.290 0.772
## time -0.001423 0.006021 27.525849 -0.236 0.815
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cnc..2
## cncl.pbvn.2 -0.368
## time -0.599 -0.006
x3<- aic(model.pubevent)
model.restr.gather <- lmer(probit.case ~ restric.gather.red +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.restr.gather)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ restric.gather.red + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 94.6 133.5 -38.3 76.6 546
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1704 -0.4927 -0.0354 0.5313 6.7138
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.235730 0.48552
## time 0.001022 0.03197 -0.65
## Residual 0.047838 0.21872
## Number of obs: 555, groups: iso_code, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.302641 0.119205 70.361538 -19.317 < 2e-16 ***
## restric.gather.red1 -0.296975 0.110053 542.504290 -2.698 0.007183 **
## restric.gather.red3 -0.218769 0.080203 523.422597 -2.728 0.006592 **
## restric.gather.red4 -0.276462 0.079468 521.500973 -3.479 0.000546 ***
## time -0.000747 0.006262 27.758479 -0.119 0.905896
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) rst..1 rst..3 rst..4
## rstrc.gth.1 -0.464
## rstrc.gth.3 -0.576 0.686
## rstrc.gth.4 -0.613 0.674 0.878
## time -0.497 -0.036 -0.035 -0.031
x4 <- aic(model.restr.gather)
model.cls.trans <- lmer(probit.case ~ close.transport.red +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.cls.trans)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ close.transport.red + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 86.7 116.9 -36.3 72.7 548
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1660 -0.5231 -0.0299 0.5296 6.9448
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.247529 0.4975
## time 0.001024 0.0320 -0.65
## Residual 0.047302 0.2175
## Number of obs: 555, groups: iso_code, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.450992 0.098912 31.484667 -24.780 < 2e-16 ***
## close.transport.red1 -0.174516 0.042036 549.166296 -4.152 3.83e-05 ***
## time -0.000759 0.006262 27.902486 -0.121 0.904
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cls..1
## cls.trnsp.1 -0.252
## time -0.628 -0.025
x5 <- aic(model.cls.trans)
model.stay.home <- lmer(probit.case ~ stay.home.red +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.stay.home)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ stay.home.red + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 103.3 133.5 -44.6 89.3 548
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2604 -0.5068 -0.0481 0.5328 7.1337
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2295284 0.47909
## time 0.0009564 0.03093 -0.64
## Residual 0.0492192 0.22185
## Number of obs: 555, groups: iso_code, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.547814 0.094101 30.050805 -27.075 <2e-16 ***
## stay.home.red1 -0.013358 0.035099 546.888338 -0.381 0.704
## time -0.001325 0.006077 27.999834 -0.218 0.829
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sty..1
## stay.hm.rd1 -0.192
## time -0.629 -0.038
x6 <- aic(model.stay.home)
model.restr.move <- lmer(probit.case ~ restric.move.red +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.restr.move)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ restric.move.red + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 101.8 132.0 -43.9 87.8 548
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3178 -0.5154 -0.0484 0.5343 7.1029
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2400995 0.49000
## time 0.0009943 0.03153 -0.65
## Residual 0.0489042 0.22114
## Number of obs: 555, groups: iso_code, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.534366 0.095680 28.964016 -26.488 <2e-16 ***
## restric.move.red2 -0.046292 0.036008 547.304999 -1.286 0.199
## time -0.001214 0.006183 27.635482 -0.196 0.846
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) rst..2
## rstrc.mv.r2 -0.165
## time -0.643 -0.025
x7 <- aic(model.restr.move)
model.inter.trav <- lmer(probit.case ~ inter.travel +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.inter.trav)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ inter.travel + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 104.3 143.2 -43.2 86.3 546
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3336 -0.4979 -0.0561 0.5121 7.1584
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2363377 0.48615
## time 0.0009839 0.03137 -0.65
## Residual 0.0488493 0.22102
## Number of obs: 555, groups: iso_code, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.474062 0.109431 49.579318 -22.608 <2e-16 ***
## inter.travel2 -0.106001 0.061352 547.033631 -1.728 0.0846 .
## inter.travel3 -0.072392 0.068242 553.914190 -1.061 0.2892
## inter.travel4 -0.088579 0.077671 550.851424 -1.140 0.2546
## time -0.001490 0.006165 27.593026 -0.242 0.8109
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) intr.2 intr.3 intr.4
## inter.trvl2 -0.465
## inter.trvl3 -0.475 0.742
## inter.trvl4 -0.457 0.717 0.696
## time -0.571 0.006 -0.003 0.044
x8 <- aic(model.inter.trav)
model.test <- lmer(probit.case ~ test.policy.red +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.test)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ test.policy.red + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 102.2 132.4 -44.1 88.2 546
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2927 -0.5068 -0.0320 0.5216 7.1756
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2306891 0.48030
## time 0.0009633 0.03104 -0.64
## Residual 0.0490835 0.22155
## Number of obs: 553, groups: iso_code, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.562544 0.096770 32.166052 -26.481 <2e-16 ***
## test.policy.red1 0.014886 0.050557 525.740203 0.294 0.769
## time -0.001444 0.006101 27.889687 -0.237 0.815
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) tst..1
## tst.plcy.r1 -0.291
## time -0.608 -0.055
x9 <- aic(model.test)
model.tracing <- lmer(probit.case ~ contact.tracing +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.tracing)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ contact.tracing + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 105.0 139.4 -44.5 89.0 542
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3282 -0.5107 -0.0292 0.5465 7.2002
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2285061 0.47802
## time 0.0009832 0.03136 -0.64
## Residual 0.0490618 0.22150
## Number of obs: 550, groups: iso_code, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.647e+00 1.439e-01 1.302e+02 -18.388 <2e-16 ***
## contact.tracing1 7.209e-02 1.133e-01 5.272e+02 0.636 0.525
## contact.tracing2 1.079e-01 1.159e-01 5.186e+02 0.931 0.352
## time -8.887e-04 6.165e-03 2.810e+01 -0.144 0.886
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cntc.1 cntc.2
## cntct.trcn1 -0.744
## cntct.trcn2 -0.762 0.930
## time -0.458 0.047 0.062
x10 <- aic(model.tracing)
model.mask <- lmer(probit.case ~ wear.mask +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.mask)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ wear.mask + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 109.1 148.0 -45.6 91.1 544
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3215 -0.5246 -0.0494 0.5305 7.1577
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2281657 0.47767
## time 0.0009572 0.03094 -0.64
## Residual 0.0494212 0.22231
## Number of obs: 553, groups: iso_code, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.621808 0.374853 29.197539 -6.994 1.05e-07 ***
## wear.mask2 0.073485 0.394960 33.419385 0.186 0.854
## wear.mask3 0.081558 0.377430 27.963831 0.216 0.830
## wear.mask4 0.063663 0.377380 27.949729 0.169 0.867
## time -0.001686 0.006118 28.510302 -0.276 0.785
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) wr.ms2 wr.ms3 wr.ms4
## wear.mask2 -0.927
## wear.mask3 -0.966 0.949
## wear.mask4 -0.968 0.952 0.993
## time -0.161 0.017 -0.007 0.006
x11 <- aic(model.mask)
model.info <- lmer(probit.case ~ info.campaign +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.info)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ info.campaign + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 96.6 131.2 -40.3 80.6 547
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0543 -0.5110 -0.0343 0.5388 7.0399
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.234285 0.48403
## time 0.001039 0.03224 -0.69
## Residual 0.048432 0.22007
## Number of obs: 555, groups: iso_code, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.665876 0.251606 465.932475 -10.595 <2e-16 ***
## info.campaign1 0.387860 0.251629 525.319857 1.541 0.124
## info.campaign2 0.105837 0.233361 506.472699 0.454 0.650
## time -0.002833 0.006328 27.825543 -0.448 0.658
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) inf.c1 inf.c2
## info.cmpgn1 -0.863
## info.cmpgn2 -0.929 0.927
## time -0.270 -0.015 0.015
x12 <- aic(model.info)
#population: no
model.pop <- lmer(probit.case ~ pop_log +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.pop)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ pop_log + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 103.3 133.5 -44.6 89.3 548
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2884 -0.5095 -0.0432 0.5308 7.1668
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2202337 0.46929
## time 0.0009495 0.03081 -0.62
## Residual 0.0492504 0.22192
## Number of obs: 555, groups: iso_code, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.264569 0.679272 28.403577 -3.334 0.00239 **
## pop_log -0.017000 0.039448 28.004982 -0.431 0.66980
## time -0.001411 0.006052 28.044731 -0.233 0.81741
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) pop_lg
## pop_log -0.991
## time -0.084 0.000
y1 <- aic(model.pop)
#population density
model.popden <- lmer(probit.case ~ pop_denlog +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.popden)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ pop_denlog + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 97.6 127.8 -41.8 83.6 548
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2930 -0.5117 -0.0337 0.5312 7.1691
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.1776080 0.42144
## time 0.0009495 0.03081 -0.62
## Residual 0.0492504 0.22192
## Number of obs: 555, groups: iso_code, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.008411 0.223225 30.842406 -8.997 3.93e-10 ***
## pop_denlog -0.107027 0.040693 27.997629 -2.630 0.0137 *
## time -0.001406 0.006053 28.044514 -0.232 0.8180
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) pp_dnl
## pop_denlog -0.931
## time -0.231 0.000
y2 <- aic(model.popden)
# median age
model.age <- lmer(probit.case ~ median_age +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ median_age + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 102.2 132.5 -44.1 88.2 548
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2915 -0.5070 -0.0349 0.5270 7.1630
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2162986 0.46508
## time 0.0009495 0.03081 -0.63
## Residual 0.0492502 0.22192
## Number of obs: 555, groups: iso_code, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.201318 0.333172 29.660521 -6.607 2.74e-07 ***
## median_age -0.011469 0.010413 28.012157 -1.101 0.280
## time -0.001411 0.006052 28.044495 -0.233 0.817
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) medn_g
## median_age -0.963
## time -0.173 0.000
y3 <- aic(model.age)
# gdp
model.gdp <- lmer(probit.case ~ gdp.log +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.gdp)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ gdp.log + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 100.9 131.2 -43.5 86.9 548
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2935 -0.5043 -0.0357 0.5266 7.1522
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2319973 0.48166
## time 0.0009494 0.03081 -0.68
## Residual 0.0492505 0.22192
## Number of obs: 555, groups: iso_code, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.526549 0.631479 28.576159 -2.417 0.0222 *
## gdp.log -0.106908 0.064952 27.999482 -1.646 0.1110
## time -0.001418 0.006052 28.045002 -0.234 0.8164
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) gdp.lg
## gdp.log -0.989
## time -0.102 0.000
y4 <- aic(model.gdp)
# test per 1000 => too many missing because the early stage of pandemic
model.test1000 <- lmer(probit.case ~ test.thousand +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.test1000)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ test.thousand + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 102.3 132.6 -44.2 88.3 548
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2922 -0.5113 -0.0432 0.5331 7.1611
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.233387 0.48310
## time 0.000932 0.03053 -0.66
## Residual 0.049246 0.22191
## Number of obs: 555, groups: iso_code, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.525e+00 9.730e-02 3.182e+01 -25.949 <2e-16 ***
## test.thousand -9.077e-05 8.569e-05 2.966e+01 -1.059 0.298
## time -6.059e-05 6.135e-03 3.011e+01 -0.010 0.992
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) tst.th
## test.thosnd -0.290
## time -0.561 -0.208
y5 <- aic(model.test1000)
# percentage of urban
model.per_urban <- lmer(probit.case ~ per_urban +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.per_urban)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ per_urban + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 100.0 130.2 -43.0 86.0 548
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2963 -0.5018 -0.0356 0.5263 7.1450
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2407254 0.49064
## time 0.0009493 0.03081 -0.71
## Residual 0.0492508 0.22193
## Number of obs: 555, groups: iso_code, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.254156 0.177031 35.178733 -12.733 9.86e-15 ***
## per_urban -0.004980 0.002481 28.003151 -2.007 0.0545 .
## time -0.001422 0.006052 28.044844 -0.235 0.8159
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) pr_rbn
## per_urban -0.846
## time -0.382 0.000
y6 <- aic(model.per_urban)
# uhc_index
model.uhc <- lmer(probit.case ~ uhc_index +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.uhc)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ uhc_index + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 103.0 133.2 -44.5 89.0 548
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2921 -0.5068 -0.0358 0.5279 7.1601
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2253242 0.47468
## time 0.0009495 0.03081 -0.64
## Residual 0.0492503 0.22192
## Number of obs: 555, groups: iso_code, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.297240 0.408576 29.177211 -5.623 4.41e-06 ***
## uhc_index -0.003871 0.005987 28.004973 -0.647 0.523
## time -0.001413 0.006052 28.044511 -0.233 0.817
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) uhc_nd
## uhc_index -0.975
## time -0.145 0.000
y7 <- aic(model.uhc)
# health worker density
model.hworker <- lmer(probit.case ~ healthworker_den +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.hworker)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ healthworker_den + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 102.0 132.2 -44.0 88.0 548
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2972 -0.5022 -0.0356 0.5291 7.1524
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2270257 0.47647
## time 0.0009495 0.03081 -0.66
## Residual 0.0492503 0.22192
## Number of obs: 555, groups: iso_code, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.404549 0.153802 35.176866 -15.634 <2e-16 ***
## healthworker_den -0.002500 0.002054 28.006164 -1.217 0.234
## time -0.001415 0.006052 28.044592 -0.234 0.817
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) hlthw_
## hlthwrkr_dn -0.802
## time -0.398 0.000
y8 <- aic(model.hworker)
#Income group
model.income <- lmer(probit.case ~ incomegroup2 +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.income)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ incomegroup2 + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 102.0 136.6 -43.0 86.0 547
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2934 -0.5038 -0.0287 0.5261 7.1545
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2221627 0.47134
## time 0.0009495 0.03081 -0.68
## Residual 0.0492501 0.22192
## Number of obs: 555, groups: iso_code, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.443163 0.115400 35.758920 -21.171 <2e-16 ***
## incomegroup2upper_income -0.102084 0.164078 28.022348 -0.622 0.539
## incomegroup2high_income -0.300860 0.157239 27.999869 -1.913 0.066 .
## time -0.001415 0.006052 28.044568 -0.234 0.817
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) incmgrp2p_ incmgrp2h_
## incmgrp2pp_ -0.498
## incmgrp2hg_ -0.519 0.365
## time -0.541 0.001 0.000
y9 <- aic(model.income)
# SDI
model.sdi <- lmer(probit.case ~ sdi_index +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.sdi)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ sdi_index + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 101.7 131.9 -43.8 87.7 548
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2950 -0.5065 -0.0335 0.5297 7.1547
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2254544 0.47482
## time 0.0009495 0.03081 -0.66
## Residual 0.0492504 0.22192
## Number of obs: 555, groups: iso_code, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.140177 0.324088 29.971522 -6.604 2.63e-07 ***
## sdi_index -0.626891 0.470193 28.002484 -1.333 0.193
## time -0.001416 0.006052 28.044655 -0.234 0.817
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sd_ndx
## sdi_index -0.959
## time -0.189 0.000
y10 <- aic(model.sdi)
# mobi.composite
model.mobi <- lmer(probit.case ~ mobi.composite +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.mobi)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ mobi.composite + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## -39.8 -10.9 26.9 -53.8 451
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2721 -0.5685 -0.0579 0.4959 3.5716
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.140904 0.37537
## time 0.001008 0.03175 -0.57
## Residual 0.036810 0.19186
## Number of obs: 458, groups: iso_code, 23
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.697e+00 8.170e-02 2.471e+01 -33.011 <2e-16 ***
## mobi.composite 7.427e-04 1.003e-03 4.491e+02 0.741 0.459
## time 5.051e-03 6.813e-03 2.179e+01 0.741 0.466
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) mb.cmp
## mobi.compst 0.193
## time -0.562 0.056
y11 <- aic(model.mobi)
# pp.vac.100
vac.100 <- lmer(probit.case ~ pp.vac.100 +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(vac.100)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ pp.vac.100 + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 102.4 145.6 -41.2 82.4 545
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2500 -0.4979 -0.0390 0.5044 7.0300
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2418162 0.49175
## time 0.0009176 0.03029 -0.65
## Residual 0.0485558 0.22035
## Number of obs: 555, groups: iso_code, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.531025 0.096004 29.241183 -26.364 < 2e-16 ***
## pp.vac.1001 -0.060768 0.040118 529.889865 -1.515 0.13043
## pp.vac.1002 -0.117374 0.058890 549.454377 -1.993 0.04675 *
## pp.vac.1003 -0.163785 0.074183 553.142389 -2.208 0.02766 *
## pp.vac.1004 -0.256442 0.097664 506.277510 -2.626 0.00891 **
## time 0.006671 0.006734 44.127606 0.991 0.32725
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) p..1001 p..1002 p..1003 p..1004
## pp.vac.1001 -0.161
## pp.vac.1002 -0.135 0.724
## pp.vac.1003 -0.108 0.684 0.806
## pp.vac.1004 -0.087 0.556 0.656 0.745
## time -0.520 -0.313 -0.383 -0.442 -0.429
y11 <- aic(vac.100)
# ful.vac.100
ful.vac.100 <- lmer(probit.case ~ ful.vac.100 +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(ful.vac.100)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ ful.vac.100 + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 91.4 121.7 -38.7 77.4 548
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3804 -0.4942 -0.0434 0.5060 7.2444
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2203538 0.46942
## time 0.0007844 0.02801 -0.60
## Residual 0.0484613 0.22014
## Number of obs: 555, groups: iso_code, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.576378 0.090747 28.201276 -28.391 < 2e-16 ***
## ful.vac.100 -0.007512 0.002124 398.790347 -3.537 0.000452 ***
## time 0.004567 0.005792 32.010772 0.788 0.436236
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) f..100
## ful.vac.100 0.068
## time -0.604 -0.292
y12 <- aic(ful.vac.100)
aic_table <- rbind(x1, x4, x5,x9 , y2, y6, y11, y12)
model.mul <- lmer(probit.case ~ close.transport.red + pp.vac.100 + school.close.red + pop_denlog + restric.gather.red +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.mul)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ close.transport.red + pp.vac.100 + school.close.red +
## pop_denlog + restric.gather.red + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 49.6 123.0 -7.8 15.6 538
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7098 -0.5495 0.0116 0.4771 6.9168
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2146046 0.46325
## time 0.0009648 0.03106 -0.64
## Residual 0.0426146 0.20643
## Number of obs: 555, groups: iso_code, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.801946 0.252927 36.922632 -7.124 1.96e-08 ***
## close.transport.red1 -0.225375 0.049499 528.896720 -4.553 6.57e-06 ***
## pp.vac.1001 -0.090380 0.037939 527.396458 -2.382 0.01756 *
## pp.vac.1002 -0.154076 0.055698 547.133298 -2.766 0.00586 **
## pp.vac.1003 -0.169693 0.070545 552.678279 -2.405 0.01648 *
## pp.vac.1004 -0.291346 0.093020 525.947776 -3.132 0.00183 **
## school.close.red2 0.013337 0.041507 538.239817 0.321 0.74810
## school.close.red3 0.184678 0.035316 535.362646 5.229 2.44e-07 ***
## pop_denlog -0.102345 0.044287 28.754591 -2.311 0.02821 *
## restric.gather.red1 -0.168341 0.118995 543.482660 -1.415 0.15773
## restric.gather.red3 -0.089218 0.091998 548.902903 -0.970 0.33258
## restric.gather.red4 -0.195709 0.089499 540.881089 -2.187 0.02919 *
## time 0.009285 0.006778 41.125757 1.370 0.17813
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
performance(model.mul)
## # Indices of model performance
##
## AIC | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma
## ------------------------------------------------------------------
## 49.611 | 123.033 | 0.828 | 0.185 | 0.789 | 0.196 | 0.206
# #====================
# model.mul <- lmer(probit.case ~ close.transport.red + pp.vac.1 + school.close.red + pop_denlog + restric.gather.red + test.policy +
# time + (time|iso_code),
# data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
# summary(model.mul)
# performance(model.mul1)
#
#
#
# anova(model.mul, model.mul1)
# # => No different select the simplifier model
#devtools::install_github("goodekat/redres")
library(redres)
data0.anal1 <- data0.anal %>% filter( is.na(pp.vac.100) ==F) %>% filter(is.na(close.transport.red) == F) %>% filter(is.na(school.close.red) == F) %>% filter(is.na(restric.gather.red) == F) %>%
mutate(pop_denlog = log(pop_density))
rc_resids <- compute_redres(model.mul)
pm_resids <- compute_redres(model.mul, type = "pearson_mar")
## Warning in (subset(re_df, re_df$grp != "Residual")$sdcor)^2 * diag(q): longer
## object length is not a multiple of shorter object length
sc_resids <- compute_redres(model.mul, type = "std_cond")
## Warning in (subset(re_df, re_df$grp != "Residual")$sdcor)^2 * diag(q): longer
## object length is not a multiple of shorter object length
resids <- data.frame(data0.anal1$iso_code, rc_resids, pm_resids, sc_resids)
plot_redres(model.mul, type = "std_cond")
## Warning in (subset(re_df, re_df$grp != "Residual")$sdcor)^2 * diag(q): longer
## object length is not a multiple of shorter object length
plot_resqq(model.mul)
plot_ranef(model.mul)
plot(data0.anal1$time, rc_resids,ylim=c(-3,3))
abline(h=0, col="blue")
plot(data0.anal1$iso_code, sc_resids)
abline(h=c(0,2.96,-2.96), col="blue")
# # run other diagnosis
# p_load(patchwork, randomForest)
#
#check_model(model.mul)
#
plot(check_distribution(model.mul))
## Warning: Removed 6 rows containing missing values (geom_segment).
## Warning: Removed 6 rows containing missing values (geom_point).
#
check_collinearity(model.mul)
## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF Increased SE Tolerance
## close.transport.red 1.54 1.24 0.65
## pp.vac.100 1.34 1.16 0.74
## school.close.red 1.22 1.10 0.82
## pop_denlog 1.02 1.01 0.98
## restric.gather.red 1.56 1.25 0.64
## time 1.25 1.12 0.80
#
# pp_check(model.mul, 250)
library(glmmTMB)
## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.3.4
## Current Matrix version is 1.3.3
## Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package
##
## Attaching package: 'glmmTMB'
## The following object is masked from 'package:VGAM':
##
## betabinomial
# model.glm <- glmmTMB(week.case.eff ~ close.transport.red + pp.vac.100 + school.close.red + pop_denlog + restric.gather.red + test.policy +
# time + (time|iso_code), data = data0.anal, family = beta_family(link="probit"), REML = TRUE)
#
# summary(model.glm)
model.glm <- glmmTMB(week.case.eff ~ close.transport.red + pp.vac.100 + school.close.red + pop_denlog + restric.gather.red +
time + (time|iso_code), data = data0.anal, family = beta_family(link="probit"), REML = TRUE)
summary(model.glm)
## Family: beta ( probit )
## Formula:
## week.case.eff ~ close.transport.red + pp.vac.100 + school.close.red +
## pop_denlog + restric.gather.red + time + (time | iso_code)
## Data: data0.anal
##
## AIC BIC logLik deviance df.resid
## -4328.8 -4255.4 2181.4 -4362.8 551
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.1080878 0.32877
## time 0.0005248 0.02291 -0.60
## Number of obs: 555, groups: iso_code, 28
##
## Dispersion parameter for beta family (): 234
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.718147 0.200878 -8.553 < 2e-16 ***
## close.transport.red1 -0.150334 0.053208 -2.825 0.004722 **
## pp.vac.1001 -0.040151 0.039477 -1.017 0.309119
## pp.vac.1002 -0.116616 0.054956 -2.122 0.033837 *
## pp.vac.1003 -0.111118 0.068877 -1.613 0.106682
## pp.vac.1004 -0.242739 0.095917 -2.531 0.011383 *
## school.close.red2 0.001774 0.043782 0.041 0.967688
## school.close.red3 0.128579 0.036852 3.489 0.000485 ***
## pop_denlog -0.084909 0.035125 -2.417 0.015635 *
## restric.gather.red1 -0.223756 0.125291 -1.786 0.074117 .
## restric.gather.red3 -0.189786 0.098053 -1.936 0.052924 .
## restric.gather.red4 -0.235337 0.092128 -2.554 0.010636 *
## time 0.006673 0.005550 1.202 0.229177
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
p_load(plotly, modelbased)
data0.anal1$Predicted <- estimate_response(model.glm)$Predicted
## Warning in get_predicted.glmmTMB(model, data = data, predict = predict, : `predict = 'prediction'` is currently not available for glmmTMB models.
## Changing to `predict = 'expectation'`.
plot2 <- data0.anal1 %>%
ggplot() +
geom_line(aes(x = log(week.case.eff), y = log(week.case.eff)), linetype = "dashed") +
geom_point(aes(x = log(week.case.eff) , y = log(Predicted), key=iso_code), color = "red") +
ylab("wADGR (predicted)") + xlab("wADGR (observed)")
## Warning: Ignoring unknown aesthetics: key
plot2
ggplotly(plot2, source = "select", tooltip = c("key") )
# Develop function marginal
amex<- function(ame){
x1 <- data.frame(list("AME" = ame*100, "factors" = as.character(substitute(ame))))
}
summary(model.glm)
## Family: beta ( probit )
## Formula:
## week.case.eff ~ close.transport.red + pp.vac.100 + school.close.red +
## pop_denlog + restric.gather.red + time + (time | iso_code)
## Data: data0.anal
##
## AIC BIC logLik deviance df.resid
## -4328.8 -4255.4 2181.4 -4362.8 551
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.1080878 0.32877
## time 0.0005248 0.02291 -0.60
## Number of obs: 555, groups: iso_code, 28
##
## Dispersion parameter for beta family (): 234
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.718147 0.200878 -8.553 < 2e-16 ***
## close.transport.red1 -0.150334 0.053208 -2.825 0.004722 **
## pp.vac.1001 -0.040151 0.039477 -1.017 0.309119
## pp.vac.1002 -0.116616 0.054956 -2.122 0.033837 *
## pp.vac.1003 -0.111118 0.068877 -1.613 0.106682
## pp.vac.1004 -0.242739 0.095917 -2.531 0.011383 *
## school.close.red2 0.001774 0.043782 0.041 0.967688
## school.close.red3 0.128579 0.036852 3.489 0.000485 ***
## pop_denlog -0.084909 0.035125 -2.417 0.015635 *
## restric.gather.red1 -0.223756 0.125291 -1.786 0.074117 .
## restric.gather.red3 -0.189786 0.098053 -1.936 0.052924 .
## restric.gather.red4 -0.235337 0.092128 -2.554 0.010636 *
## time 0.006673 0.005550 1.202 0.229177
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coef_model <- summary(model.glm)$coefficients$cond[,1][-1]
coef <- as.data.frame(coef_model)
# AME time
AME.time <- coef_model[12] * mean(dnorm(predict(model.glm, type="link")))
z12 <- amex(AME.time)
# pop.density
AME.popden<- coef_model[8]* mean(dnorm(predict(model.glm, type="link")))
z8 <- amex(AME.popden)
# vac
vac0 <- predict(model.glm, newdata=data.frame(list(
"close.transport.red"=data0.anal1$close.transport.red,
"pp.vac.100"="0",
"school.close.red"=data0.anal1$school.close.red,
"pop_denlog"=data0.anal1$pop_denlog,
"restric.gather.red"= data0.anal1$restric.gather.red,
"test.policy"=data0.anal1$test.policy,
"time"=data0.anal1$time,
"iso_code"=data0.anal1$iso_code) ), type = "response")
vac1 <- predict(model.glm, newdata=data.frame(list(
"close.transport.red"=data0.anal1$close.transport.red,
"pp.vac.100"="1",
"school.close.red"=data0.anal1$school.close.red,
"pop_denlog"=data0.anal1$pop_denlog,
"restric.gather.red"= data0.anal1$restric.gather.red,
"test.policy"=data0.anal1$test.policy,
"time"=data0.anal1$time,
"iso_code"=data0.anal1$iso_code) ), type = "response")
vac2 <- predict(model.glm, newdata=data.frame(list(
"close.transport.red"=data0.anal1$close.transport.red,
"pp.vac.100"="2",
"school.close.red"=data0.anal1$school.close.red,
"pop_denlog"=data0.anal1$pop_denlog,
"restric.gather.red"= data0.anal1$restric.gather.red,
"test.policy"=data0.anal1$test.policy,
"time"=data0.anal1$time,
"iso_code"=data0.anal1$iso_code) ), type = "response")
vac3 <- predict(model.glm, newdata=data.frame(list(
"close.transport.red"=data0.anal1$close.transport.red,
"pp.vac.100"="3",
"school.close.red"=data0.anal1$school.close.red,
"pop_denlog"=data0.anal1$pop_denlog,
"restric.gather.red"= data0.anal1$restric.gather.red,
"test.policy"=data0.anal1$test.policy,
"time"=data0.anal1$time,
"iso_code"=data0.anal1$iso_code) ), type = "response")
vac4 <- predict(model.glm, newdata=data.frame(list(
"close.transport.red"=data0.anal1$close.transport.red,
"pp.vac.100"="4",
"school.close.red"=data0.anal1$school.close.red,
"pop_denlog"=data0.anal1$pop_denlog,
"restric.gather.red"= data0.anal1$restric.gather.red,
"test.policy"=data0.anal1$test.policy,
"time"=data0.anal1$time,
"iso_code"=data0.anal1$iso_code) ), type = "response")
vac10<- vac1 - vac0
vac20<-vac2-vac0
vac30<-vac3-vac0
vac40<-vac4-vac0
AME.vac10 <- mean(vac10)
AME.vac20<-mean(vac20)
AME.vac30<-mean(vac30)
AME.vac40<-mean(vac40)
z2<-amex(AME.vac10)
z3<-amex(AME.vac20)
z4<-amex(AME.vac30)
z5<-amex(AME.vac40)
#
# AME.vac <- coef_model[2]* mean(dnorm(predict(model.glm, type="link")))
# z2 <- amex(AME.vac)
# transport
trs0 <- predict(model.glm, newdata=data.frame(list(
"close.transport.red"="0",
"pp.vac.100"=data0.anal1$pp.vac.100,
"school.close.red"=data0.anal1$school.close.red,
"pop_denlog"=data0.anal1$pop_denlog,
"restric.gather.red"= data0.anal1$restric.gather.red,
"test.policy"=data0.anal1$test.policy,
"time"=data0.anal1$time,
"iso_code"=data0.anal1$iso_code) ), type = "response")
trs1 <- predict(model.glm, newdata=data.frame(list(
"close.transport.red"="1",
"pp.vac.100"=data0.anal1$pp.vac.100,
"school.close.red"=data0.anal1$school.close.red,
"pop_denlog"=data0.anal1$pop_denlog,
"restric.gather.red"= data0.anal1$restric.gather.red,
"test.policy"=data0.anal1$test.policy,
"time"=data0.anal1$time,
"iso_code"=data0.anal1$iso_code) ), type = "response")
trs10 <- trs1 - trs0
AME.trs10 <- mean(trs10)
z1 <- amex(AME.trs10)
#school close
sch0 <- predict(model.glm, newdata=data.frame(list(
"close.transport.red"=data0.anal1$close.transport.red,
"pp.vac.100"=data0.anal1$pp.vac.100,
"school.close.red"="0",
"pop_denlog"=data0.anal1$pop_denlog,
"restric.gather.red"= data0.anal1$restric.gather.red,
"test.policy"=data0.anal1$test.policy,
"time"=data0.anal1$time,
"iso_code"=data0.anal1$iso_code) ), type = "response")
sch2 <- predict(model.glm, newdata=data.frame(list(
"close.transport.red"=data0.anal1$close.transport.red,
"pp.vac.100"=data0.anal1$pp.vac.100,
"school.close.red"="2",
"pop_denlog"=data0.anal1$pop_denlog,
"restric.gather.red"= data0.anal1$restric.gather.red,
"test.policy"=data0.anal1$test.policy,
"time"=data0.anal1$time,
"iso_code"=data0.anal1$iso_code) ), type = "response")
sch3 <- predict(model.glm, newdata=data.frame(list(
"close.transport.red"=data0.anal1$close.transport.red,
"pp.vac.100"=data0.anal1$pp.vac.100,
"school.close.red"="3",
"pop_denlog"=data0.anal1$pop_denlog,
"restric.gather.red"= data0.anal1$restric.gather.red,
"test.policy"=data0.anal1$test.policy,
"time"=data0.anal1$time,
"iso_code"=data0.anal1$iso_code) ), type = "response")
sch20 <- sch2 - sch0
sch30 <- sch3 - sch0
AME.sch20 <- mean(sch20)
AME.sch30 <- mean(sch30)
z6 <- amex(AME.sch20)
z7 <- amex(AME.sch30)
# restrict gather
gather0 <- predict(model.glm, newdata=data.frame(list(
"close.transport.red"=data0.anal1$close.transport.red,
"pp.vac.100"=data0.anal1$pp.vac.100,
"school.close.red"=data0.anal1$school.close.red,
"pop_denlog"=data0.anal1$pop_denlog,
"restric.gather.red"= "0",
"test.policy"=data0.anal1$test.policy,
"time"=data0.anal1$time,
"iso_code"=data0.anal1$iso_code) ), type = "response")
gather1 <- predict(model.glm, newdata=data.frame(list(
"close.transport.red"=data0.anal1$close.transport.red,
"pp.vac.100"=data0.anal1$pp.vac.100,
"school.close.red"=data0.anal1$school.close.red,
"pop_denlog"=data0.anal1$pop_denlog,
"restric.gather.red"= "1",
"test.policy"=data0.anal1$test.policy,
"time"=data0.anal1$time,
"iso_code"=data0.anal1$iso_code) ), type = "response")
gather3 <- predict(model.glm, newdata=data.frame(list(
"close.transport.red"=data0.anal1$close.transport.red,
"pp.vac.100"=data0.anal1$pp.vac.100,
"school.close.red"=data0.anal1$school.close.red,
"pop_denlog"=data0.anal1$pop_denlog,
"restric.gather.red"= "3",
"test.policy"=data0.anal1$test.policy,
"time"=data0.anal1$time,
"iso_code"=data0.anal1$iso_code) ), type = "response")
gather4 <- predict(model.glm, newdata=data.frame(list(
"close.transport.red"=data0.anal1$close.transport.red,
"pp.vac.100"=data0.anal1$pp.vac.100,
"school.close.red"=data0.anal1$school.close.red,
"pop_denlog"=data0.anal1$pop_denlog,
"restric.gather.red"= "4",
"test.policy"=data0.anal1$test.policy,
"time"=data0.anal1$time,
"iso_code"=data0.anal1$iso_code) ), type = "response")
gather10 <- gather1 - gather0
gather30 <- gather3 - gather0
gather40 <- gather4 - gather0
AME.gather10 <- mean(gather10)
AME.gather30 <- mean(gather30)
AME.gather40 <- mean(gather40)
z9 <- amex(AME.gather10)
z10 <- amex(AME.gather30)
z11 <- amex(AME.gather40)
# Test policy
# test0 <- predict(model.glm, newdata=data.frame(list(
# "close.transport.red"=data0.anal1$close.transport.red,
# "pp.vac.100"=data0.anal1$pp.vac.100,
# "school.close.red"=data0.anal1$school.close.red,
# "pop_denlog"=data0.anal1$pop_denlog,
# "restric.gather.red"= data0.anal1$restric.gather.red,
# "test.policy"="0",
# "time"=data0.anal1$time,
# "iso_code"=data0.anal1$iso_code) ), type = "response")
# test1 <- predict(model.glm, newdata=data.frame(list(
# "close.transport.red"=data0.anal1$close.transport.red,
# "pp.vac.100"=data0.anal1$pp.vac.100,
# "school.close.red"=data0.anal1$school.close.red,
# "pop_denlog"=data0.anal1$pop_denlog,
# "restric.gather.red"= data0.anal1$restric.gather.red,
# "test.policy"="1",
# "time"=data0.anal1$time,
# "iso_code"=data0.anal1$iso_code) ), type = "response")
# test2 <- predict(model.glm, newdata=data.frame(list(
# "close.transport.red"=data0.anal1$close.transport.red,
# "pp.vac.100"=data0.anal1$pp.vac.100,
# "school.close.red"=data0.anal1$school.close.red,
# "pop_denlog"=data0.anal1$pop_denlog,
# "restric.gather.red"= data0.anal1$restric.gather.red,
# "test.policy"="2",
# "time"=data0.anal1$time,
# "iso_code"=data0.anal1$iso_code) ), type = "response")
# test3 <- predict(model.glm, newdata=data.frame(list(
# "close.transport.red"=data0.anal1$close.transport.red,
# "pp.vac.100"=data0.anal1$pp.vac.100,
# "school.close.red"=data0.anal1$school.close.red,
# "pop_denlog"=data0.anal1$pop_denlog,
# "restric.gather.red"= data0.anal1$restric.gather.red,
# "test.policy"="3",
# "time"=data0.anal1$time,
# "iso_code"=data0.anal1$iso_code) ), type = "response")
#
# test10 <- test1- test0
# test20 <- test2-test0
# test30 <- test3-test0
#
# AME.test10 <- mean(test10)
# AME.test20 <-mean(test20)
# AME.test30 <-mean(test30)
#
# z12 <- amex(AME.test10)
# z13 <- amex(AME.test20)
# z14 <- amex(AME.test30)
z_merge <- rbind(z1,z2,z3,z4,z5,z6,z7,z8,z9,z10,z11,z12)
term3 <- c(
"Closing public transport",
"1-<5%",
"5-<10%",
"10-<30%",
"Vaccine coverage\n>=30%",
"School close at some levels",
"School close at all levels",
"Log of population density",
"Restrictions on gathering\n100-<1000",
"Restrictions on gathering\n10-<100",
"Restrictions on gathering\n<10",
"Time\n"
)
p_load(broom.mixed)
sum_table <- tidy(model.mul, conf.int = T)
# create new table summary of model
new_table <- sum_table[2:13,]
new_table1 <- cbind(new_table, term3, z_merge)
new_table2 <- new_table1 %>% select(term, term3, AME, estimate, conf.low, conf.high)
new_table2$term3 <- factor(new_table2$term3, levels = new_table2$term3)
p.vac.period<- ggplot(new_table2, aes(x = term3, y = estimate,
ymin = conf.low, ymax = conf.high)) +
geom_hline( yintercept = 0, color = 'Black' ) + ylim(-0.6,0.6) +
geom_linerange(size =1, color = 'steelblue') + geom_point(size=1.5, color ='steelblue') + coord_flip() +
ylab("Probit Transformation of Average Weekly Growth Rate ") +
xlab("Predictors") +
theme_minimal() +
theme(strip.background = element_rect(NA),
axis.text = element_text(size = 20),
axis.title = element_text(size = 20),
legend.text = element_text(size = 10),
legend.title = element_text(size = 10),
strip.text = element_text(size = 8),
legend.position = "top")
p.vac.period
ame.plot <- ggplot(data = new_table2, aes(x= term3, y = AME)) + geom_bar(stat="identity",fill="steelblue") + coord_flip() +
geom_hline( yintercept = 0, color = 'Black' ) +
geom_text(aes(label=round(AME,2)), hjust= 1.05, color="black", size=5)+
ylab("Average marginal effect (%)") +
ylim(-1,1) +
xlab("") +
theme_minimal() +
theme(strip.background = element_rect(NA),
axis.text.y=element_blank(),
axis.text = element_text(size = 20),
axis.title = element_text(size = 20),
legend.text = element_text(size = 10),
legend.title = element_text(size = 10),
strip.text = element_text(size = 8),
legend.position = "top")
ame.plot
cairo_pdf(filename = "F:/STUDY HCA IN TAIWAN/0. ideas for thesis/COVID19 and NPIs/Data analysis/results/plot_sep/fig3_2021_vac_cate_Sep29.pdf",
width = 20, height = 10, #inch
onefile = TRUE, family = "Segoe UI")
merge_bar <- ggarrange(p.vac.period, ame.plot, ncol = 2, nrow = 1, widths = c(1.5,1))
merge_bar
# cairo_pdf(filename = "F:/STUDY HCA IN TAIWAN/0. ideas for thesis/COVID19 and NPIs/Data analysis/results/plot_sep/fig3_2021_vaccine_Sep29.pdf",
# width = 20, height = 10, #inch
# onefile = TRUE, family = "Segoe UI")
# merge_bar <- ggarrange(p.vac.period, ame.plot, ncol = 2, nrow = 1, widths = c(1.5,1))
# merge_bar
dev.off()
## png
## 2